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Abstract. We calculate the three bulk viscosity coefficients as arising from the collisions 
among phonons in superfluid neutron stars. We use effective field theory techniques to 
extract the allowed phonon collisional processes, written as a function of the equation of 
state of the system. The solution of the dynamical evolution of the phonon number density 
allows us to calculate the bulk viscosity coefficients as function of the phonon collisional rate 
and the phonon dispersion law, which depends on the neutron pairing gap. Our method of 
computation is rather general, and could be used for different superfluid systems, provided 
they share the same underlying symmetries. We find that the behavior with temperature of 
the bulk viscosity coefficients is dominated by the contributions coming from the collinear 
regime of the 2 -f-)- 3 phonon processes. For typical star radial pulsation frequencies of 
oj ~ 10^s~^ we see that the static and frequency-dependent bulk viscosity coefficients due to 
phonons are distinguishable for n > Atiq depending on the model used for the neutron gap. 
Compared to previous results from URCA and modified URCA reactions, we conclude that 
at T ~ lO^K phonon collisions give the leading contribution to the bulk viscosities in the 
core of the neutron stars, except for n ~ 2no when the opening of the URCA processes takes 
place. 
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1 Introduction 

There are strong theoretical reasons to believe that superfluidity occurs in both the inner 
crust and in the core of neutron stars [1, 2]. The fact that part of the nucleon-nucleon 
interaction is attractive for some values of the nucleonic density implies the formation of 
Cooper pairs of neutrons and Cooper pair of protons inside the star, leading to superfluidity 
and superconductivity, respectively. On the other hand, superfluidity seems to be needed to 
explain a variety of different neutron star phenomena, such as the existence of pulsar glitches 
[3], or the cooling of the star [4]. 

Superfluidity would also affect many other aspects of the neutron star dynamics, the 
corresponding hydrodynamics being essentially different from that of a normal fluid, hi 
particular, rotational and vibrational properties of the star, or the dynamics of the stars 
oscillations, would be different for a star with or without a superfluid core. Further, the 
dynamics and damping time scales of both the rotational, vibrational or oscillations modes 
of the star are governed by the value of the different transport coefficients of the matter 
inside the star [5]. In general, one may expect that the value of the different transport 
coefficients (viscosities, conductivities, etc) might attain much lower values in the superfluid 
rather than in the normal phase of the nuclear matter. Several computations of the main 
transport coefficients in the two phases of nuclear matter can be found in the literature (see, 
e.g., Refs. [6-13]). 

It is also important to recognize that the hydrodynamic equations governing the bulk 
fluctuations of a superfluid are essentially different from those of a normal fluid. At non- 
vanishing temperature one has to employ the two-fluid description of Landau [14] , which takes 
into account the motion of both the superfluid and of the normal component of the system. In 
order to describe the different dissipative processes one introduces more transport coefficients 
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than in a normal fluid. In particular, one has three independent bulk viscosities, C11C21C31 
as well as the shear viscosity and the thermal conductivity [14, 15]. The hydrodynamics 
describing the fluid motions inside a star are even much more complicated, because, apart 
from the superfluid hydro dynamical equations, one has to describe as well the fluid motions 
of the charged components (mainly, protons and electrons) of the star. 

In this article we compute the superfluid phonon contribution to the three different 
bulk viscosities in superfluid neutron stars. The phonon contribution to the shear viscosity 
in superfluid neutron stars was already computed in Ref. [13, 16]. The phonon is a collective 
mode which appears in all superfluids, and which is essential to explain the property of 
superfluidity, as found out by Landau. In field theoretical language one views the superfluid 
phonon as the Goldstone mode which appears due to the fact that the neutron condensate 
spontaneously breaks the U{1) baryonic symmetry. At low densities, the neutrons pair in a 
^Sq channel within the star, while at higher densities a '^^2 channel is preferred [1], but the 
superfluid phonon exists in both of these two different superfluid phases. 

The bulk viscosity coefficients depend on a collisional rate of phonon number changing 
scatterings. Effective field theory (EFT) techniques can be applied to study the phonon self- 
interactions [17, 18]. Then, it is seen that, at leading order (LO) in a derivative expansion, 
the phonon self-couplings are determined by the equation of state (EoS) of the superfluid. 
In this paper we consider a simplified model of neutron star made up by neutrons, protons 
and electrons, using a causal parametrization of the Akmal, Pandharipande and Ravenhall 
[19] (APR for short) EoS [2] to describe the /3-stable nuclear matter inside the star, and get 
from it all the phonon self-couplings. 

Starting from the LO phonon EFT one sees that, in principle, there are several possible 
phonon number changing processes. However, the energy and momentum conservation laws 
put kinematical constraints on the possible phonon collisions when corrections to the LO 
dispersion law are considered. If the correction to the linear dispersion law curves upward, 
one phonon can decay into two, but this process is kinematically forbidden in the opposite 
case. If the correction to the linear dispersion laws curves downward, then the most important 
process is dominated by 2 o 3 small angle scatterings, as we will discuss at length. The 
phonon dispersion law can be computed through a matching procedure with the underlying 
nucleonic microscopic theory. For neutrons pairing in a ^Sq channel, the phonon dispersion 
law curves downward [20] , and in this article we will assume that this is the case in the whole 
range of densities of the star. Thus we will assume that the collisions determining the values 
of the bulk viscosities are 2 -f-)- 3 phonon scatterings. 

As the superfluid phonon is a massless mode, one may expect that at sufficiently low 
T it might give the leading contribution to the bulk viscosity coefficients. We compare our 
results for (2 of those obtained from the contribution of both the direct URCA and modified 
URCA processes [8, 9], finding that this is indeed the case. The values of the coefficients Ci 
and (^3 as arising from both direct and modified URCA processes have been analyzed as well 
in Ref. [21]^. 

This paper is structured as follows. In Section 2 we present the EFT for superfluid 
phonons at leading order and comment on how corrections at next-to-leading order change 
the phonon dispersion law. We also present the nucleonic EoS used in this work, which is a 

^We note, however, that the computations of Ref. [21] assumed relativistic hydrodynamics instead of the 
non-relativistic domain that is considered in most of the computations of transport coefficients in the neutron 
star core. The first and third bulk viscosity coefficients of Ref. [21] do not have the same dimensions than 
those of the non-relativistic superfluid hydrodynamics. 
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common benchmark for all EoS in neutron star matter. In Section 3 we review the expressions 
for the phonon contribution to the (static) bulk viscosity coefficients, and get also expressions 
for the frequency dependent coefficients, while in Section 4 we show the phonon decay rate, 
which consists of 2 -H- 3 collisions, that is relevant for the computation of the bulk viscosity 
coefficients. Our results for the bulk viscosity coefficients are given in Section 5 and our 
conclusions in Section 6. Finally, in Appendix A we discuss the allowed kinematics when a 
phonon dispersion law curves downward and in Appendix B we give details regarding the 
phase space integral for the calculation of the phonon decay rate. We use of natural units 
h = c = kh = 1 in all intermediate computations, except in the plots as we report our results 
in CGS units. 

2 Superfluid phonons interactions and the Equation of State 

In this Section we first review the Lagrangian describing the phonon self-interactions using 
EFT techniques. At LO in a derivative expansion, all the phonon self-couplings can be 
parametrized in terms of the speed of sound, the density of the superfluid, and the derivatives 
of the speed of sound with respect to the density at zero temperature T = (Sec. 2.1). This 
is a result which is valid for any superfluid sharing the same global symmetries. Thus, the 
phonon physics at LO is determined by the T = EoS of the superfluid. In Sec. 2.2 we present 
the EoS that we use for the superfluid matter in neutron stars for the explicit computations 
of the bulk viscosity coefficients in the remaining part of this article. 

2.1 Effective field theory for superfluid phonons at leading order 

The superfluid phonon is the Goldstone mode associated to the spontaneous symmetry break- 
ing of a U{1) symmetry, which corresponds to particle number conservation. EFT techniques 
can be used to write down the effective Lagrangian associated to the superfluid phonon. The 
effective Lagrangian is then presented as an expansion in derivatives of the Goldstone field, 
the terms of this expansion being restricted by symmetry considerations. The coefficients of 
the Lagrangian can be in principle computed from the microscopic theory, through a stan- 
dard matching procedure, and thus they depend on the short range physics of the system 
under consideration. 

It has been known for a while that the leading-order term Lagrangian of the Goldstone 
mode in a superfluid system is entirely fixed by the EoS [22]. In recent publications [17, 18] 
it has been realized that at lowest order in a derivative expansion the Lagrangian reads [18] 



where P{fi) and /U are the pressure and chemical potential, respectively, of the superfluid at 
T = 0. The quantity (f is the phonon field and m is the mass of the particles that condense. 
After a Legendre transformation, one can associate this formulation to the one due to Popov 
[18, 22]. The associated Hamiltonian has also the same form as the one used by Landau to 
obtain the self-interactions of the phonons of "^He [15, 18]. 

The origin of this particular form for the LO Lagrangian is that the effective action 
associated to the theory at its minimum for constant classical field configurations has to be 



Clo = P{X) 



X = i^L- dtif 




2m 
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equal to the pressure [17]. This formulation turns out to be very advantageous, as it allows 
one to derive all the phonon properties at lowest order in a momentum expansion based on 
the knowledge of the zero temperature pressure of the superfluid. In particular, one can 
easily get the phonon dispersion law and the form of the leading phonon self-interactions, as 
well as their leading contribution to different physical processes. 

In order to see that, we expand the P{X) around ^, and after the field redefinition 



(2.2) 



to have the kinetic term canonically normalized, one can write [23] 



J~-LO 



+A 



+ ... 



(2.3) 



We have neglected above an irrelevant constant and a total time derivative term, which is 
only needed to study vortex configurations. 

The different phonon self-couplings of eq. (2.3) can be expressed as different ratios of 
derivatives of the pressure with respect to the chemical potential [23]. In particular, after 
using the thermodynamic relation dP = ^dfi, where p is the mass density at T = 0, the 
phonon velocity is 



Vph 



dP 
dfj. 



m 
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that is, it can be identified with the speed of sound at T = 0, as it is expected in the low 
momentum limit. The dispersion law obtained from this Lagrangian at tree level is exactly 



i,p = CsP. 

Defining the quantities 



p dcs 
Cs op 



w 



p d'^Cs 
Cs dp^ 



(2.5) 



we can obtain the three and four phonon self-couplings in terms of the speed of sound, the 
mass density and derivatives of the speed of sound with respect to the mass density [23] 



1 - 2u 



l-2u 
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1 - 2u{A - 5n) - 2wp 
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1 - 2n(4 - bu) - 2wp 
24^!p ' 

1 - 2u(4 - 5n) - 2wp 



(2.6) 



A next-to-leading order (NLO) Lagrangian in a derivative expansion can be constructed as 
well (see, for example, the expression of >Cnlo for the cold Fermi gas in the unitarity limit 
[18]). For our purposes, we will not need it, as we will simply compute the leading order T 
corrections of the bulk viscosity coefficients. It is however relevant for our discussion that at 
LO the phonon dispersion law is linear in the momentum, and it suffers corrections when one 
goes beyond the LO expansion. More particularly, at NLO the phonon dispersion law reads 



Ep = Csp{l +7P^ 



(2.7) 
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The sign of 7 determines whether the decay of one phonon into two is kinematically allowed 
or not (see Appendix A for more explicit details). Only dispersion laws that curve upward 
can allow such processes. The possibility of having these phonon decay processes is important 
for the value of the different transport coefficients of the superfluid. 

The value of 7 can be computed through a matching procedure with the underlying 
microscopic theory. For neutrons paring in a ^Sq channel within neutron stars it can be seen 
that [20] 

(2.8) 



' 45A2 ' 

with vp being the Fermi velocity and A the value of the gap in the ^So phase [1]. We will 
assume that 7 takes this same value in the phase, with A being the angular averaged 
value of the gap in that phase. Explicit values for the gap function used in this work will be 
provided in Sec. 4.1. 

Thus, considering that 7 < the first allowed phonon scattering will be binary colli- 
sions. In the present work we aim at computing the bulk viscosities of superfluid neutron 
star matter. Then, phonon number changing processes are needed and the first ones that 
contribute to the bulk viscosities are 2 -f-)- 3 phonon collisions. 

2.2 Equation of state for superfluid matter in neutron stars 

The speed of sound at T = as well as the different phonon self-couplings are determined by 
the EoS for neutron matter in neutron stars. A common benchmark for a nucleonic equation 
of state is the one obtained by APR. Later on Heiselberg and Hjorth-Jensen [2] parametrized 
the APR EoS of nuclear matter in a causal simple form, which will be subsequently used in 
this manuscript. The effect of neutron pairing is not considered because it is not expected 
to have a big impact in the EoS, being typically of order A^Z/x^, a quantity which remains 
small. 

The parametrized binding energy per nucleon (E/A) in nuclear matter reads 

E/A = £^yy—lzA + Soy^{l - 2x,f. (2.9) 
l + dy 

Here y is the ratio of the nucleon particle density (n) to nuclear saturation density (tiq = 
0.16 fm~'^), y = n/no, and Xp = np/riQ is the proton fraction. The nucleon density is given 

by 

n = v — ^ , (2.10) 







(2vr)= 



with pf being the Fermi momentum and v the degeneracy factor. In nuclear matter the 
degeneracy factor v is 4. The binding energy per nucleon at saturation density excluding 
Coulomb energies is £q = 15.8 MeV and the parameter 5 = 0.2 was determined by fitting the 
energy per nucleon at high density to the EoS of APR [19] with three-body forces and boost 
corrections, but taking the corrected values from Table 6 of [19]. For the symmetry energy 
at saturation density, Heiselberg et al. obtained 6*0 = 32 MeV and /3 = 0.6 for the best fit. 
The EoS is given by 

£{n, Xp) = {m + E/A{n, Xp)) , (2.11) 
with m being the mass of the nucleon, while the corresponding nucleonic energy density is 

eN{n,Xp) = £{n,Xp)n . (2.12) 
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For neutron star matter made of neutrons, protons and electrons, the total energy 
density is the sum of the nucleonic contribution (neutrons and protons), ejsf, and the one for 
electrons, eg, 

e(n, Xp, He) = eN{n, Xp) + ee(ne), (2.13) 
with He being the density of electrons. The pressure includes also both contributions 

P(n, Xp, Tie) = PN{n, Xp) + Peine) , (2.14) 
where the nucleonic and electronic contributions to the pressure are given by 
PN{n,Xp) = iiri{n,Xp) {1 - Xp) n + fj,p{n,Xp) XpU- eNin,Xp), 

Peine) = He{ne)ne - €.e{ne) , (2.15) 

being /ij the chemical potential of each specie. Those chemical potentials are calculated as 

^ _ deN{n,Xp) _ deN{n,Xp) 

fijiyn^ Xp) — , fip[n, Xp) 



dnn ' P ^ P ' 

//e(ne) ~ (STr^ne)^/^^. (2.16) 

where n„ is the density of neutrons, n„ = (1 — Xp)n. 

Nucleons in neutron stars are in /3-equilibrium against weak decay processes. The con- 
straints imposed by chemical equilibrium and charge neutrality for matter made of neutrons, 
protons and electrons are 

Pp = Pe ■ (2.17) 

These conditions fix the proportion of protons, neutrons and electrons for each particle density 
and, thus, the value of the chemical potential, energy density and pressure exerted for each 
specie at a given density. 

For the computation of the speed of sound and the different three and four phonon 
self-couplings, we will proceed as follows. The mass density is related to the energy density 
by the Einstein relation, e = p. Then, to compute the mass density we only take into account 
the nucleonic part, as m ^ mg. Further, in /3-equilibrated matter, n ~ n„. Thus, we will 
assume that the speed of sound can be computed as 



dPN{n,Xp) dPN{n,Xp) dn^ 

Csin,Xp) = i- — ^ r, 2.18 

y dpN{n,Xp) y dn„ dpN{n,Xp) 

and, for the phonon self-couplings in eq. (2.6), we will use the same chain rule. 

3 Phonon contribution to the bulk viscosity coefficients in a generic su- 
perfluid 

In this Section we review the expressions for the phonon contribution to the (static) bulk 
viscosity coefficients in a generic superfluid [15]. We then generalize these expressions for the 
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situation when there is a periodic perturbation in the system, such that the bulk viscosities 
depend on the frequency of the perturbation. 

The bulk viscosities enter as coefficients in the dissipative hydrodynamic equations. 
While the physical meaning of ^2 is the same as in a normal fluid, Cii Cs) C4 refer to dissipative 
processes which lead to entropy production only in the presence of a space-time dependent 
relative motion between the superfluid and normal fluid components. The friction forces due 
to bulk viscosities can be understood as drops in the main driving forces acting on the normal 
and superfluid components. These forces are given by the gradients of P and /i, respectively. 
One can can write that 



where Peq and /^eq are the pressure and chemical potential in the absence of bulk viscosities, 
Vn and Vg are the velocities of the normal and superfluid components, respectively, and ps is 
the mass density of the superfluid component. There are some fundamental restrictions on 
the values of these coefficients [15]. The Onsager symmetry principle imposes that Ci = C4) 
and positive entropy production requires that C2, Cs > and that Cf < (2(3- 

Let us first consider a superfluid system at finite but low temperature T, that is slightly 
away from thermodynamical equilibrium. One can extract the values of the bulk viscosities 
following a method developed by Khalatnikov [15] , which consists of studying the dynamical 
evolution of the phonon number density Np^^. It has been shown in Ref. [24] that this method 
is equivalent to the method of computing the bulk viscosities using a Boltzmann equation 
for the phonons in the relaxation time approximation. 

When a perturbation applied to a superfluid system determines a change of the number 
of phonons per unit volume, A'ph, coUisional processes tend to restore the equilibrium value 
of this quantity. The evolution equation for the phonon number can be written as 



where the rate of change is expressed as a power expansion in a "fake" phonon chemical po- 
tential, dfiph, and the decay rate of phonon changing processes, Fph. Expressing the phonon 
number as a function of the density and of entropy, NpY^{p, S), and using the linearized conti- 
nuity hydrodynamic equations for these quantities, one obtains the phonon chemical potential 
in terms of the different dissipative flows that appear in the hydrodynamic equations. These 
terms modify the equilibrium pressure and chemical potential, and with them, and making 
use of eqs. (3.1,3.2), one identifies the different bulk viscosity coefficients. 



For small departures from equilibrium and for small values of and v„, it turns out 
that [15] 



P = Pcq - Cldiv(ps(Vn - Vs)) - C2divVn 

— = — - C3div(ps(vn - Vs)) - C4divv, 



(3.1) 
(3.2) 




(3.3) 




i = 1,2,3,4 



(3.4) 




1 ) 



(3.5) 



(3.6) 
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The quantities /i,/2 have been computed in Ref. [24] for a generic superfluid, reahzing 
that in order to have non-vanishing values of the three coefficients of the bulk viscosities one 
needs to consider the phonon dispersion law beyond linear order. This result is in agreement 
with that one found out by Khalatnikov and Chernikova [25] for the phonons of ^He. In 
Ref. [24], and considering the NLO phonon dispersion eq. (2.7), it has been found that 



60r5 2 



^2 - -^(A(3)-7C(5)) (2i.c. + 3p(c.f -i.^)) , (3.7) 

where B = Cgj and C('^) is the Riemann zeta function. 

For astrophysical applications it is however more important to compute the bulk vis- 
cosity coefficients when the perturbation that leads the system out of equilibrium is periodic 
in time. Then one assumes a time evolution of the phonon density number of the form 



Nph = iVph + 5? (5Afphe*"* j , (3.8) 

where uj is the frequency of the perturbation, K denotes the real part, A^ph stands for the 
equilibrium value of the phonon density, and (5A^ph is the out of equilibrium fluctuation. A 
similar dependence is assumed for the remaining hydro dynamical variables. It is then easy 
to generalize the expressions for the transport coefficients in this situation (see for example, 
Ref. [26]). Then the bulk viscosity coefficients turn out to be complex functions which 
depend on the frequency of the perturbation. Only their real part contributes to the energy 
dissipation of the system, and this reads 

Qiio) = ^- ri^C^, i = 1,2,3,4. (3.9) 

From the above expressions one can define the value of a characteristic frequency Wc for the 
phonon collisions, defined as 

^c = y^^- (3.10) 

1 dn dfj, 

In the limit where a; <C Wc one recovers the static bulk viscosity coefficients of eqs. (3.4). 

In this article we will use the generic expressions showed in this Section to compute the 
phonon contribution to the (frequency) dependent bulk viscosity coefficients for superfluid 
neutron stars. In such a case, the values of p, n and /z should correspond to the values of 
the mass density, particle density and chemical potential of the superfluid neutrons. All the 
relevant quantities that enter into eqs. (3.9) can be computed with the EoS and the phonon 
dispersion law discussed in Sec. 2, except for the phonon decay rate, which we compute in 
the following Section. 



4 Phonon decay rate 

In this Section we compute the phonon decay rate relevant for the computation of the bulk 
viscosity coefficients. After an expansion or rarefaction of the superfluid, the system goes 
back to equilibrium after a change in the number of particles. Assuming that the phonons 
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are the relevant degrees of freedom in the T regime we consider, this imphes that we need to 
compute a colHsional rate of a process that changes the number of phonons. As discussed at 
length in Appendix A, for superfluids where the phonons have a dispersion law with negative 
values of 7, the first kinematically allowed scattering consists of 2 -H- 3 collisions. The decay 
rate associated to these collisions is given by 

Tph = I d<^5{Pa, Pb-, Pd, Pe, Pf)\\Aff{E,)f{Eh) (1 + f{Ed)) (1 + f{E,)) (1 + f{Ef)) , 

(4.1) 

where f{E) = (e^^^ — l) is the Bose-Einstein distribution function. The phase space is 
defined as 

d^biPa, Pb;Pd, Pe, Pf) = 



j=d,ej 



j=d,e,f 



k=a,b, 
d,e,f 



and in Appendix B we specify our particular choice of phase space variables. The scattering 
amplitude A describes 2 3 collisions. We compute this rate using the LO Lagrangian, 
Clo, given in eq. (2.3). 

The computation is rather involved and has to be done numerically, as one has to 
consider all the Feynman diagrams depicted in fig. 1 and fig. 2 to evaluate Tph. We give 
here some of the details and subtleties associated to the computation, to then present the 
numerical results obtained for Tph- 

We have classified all possible Feynman diagrams into two groups: those that are formed 
with one 3-phonon vertex and one 4~phonon vertex of Clo, that we call type I diagrams 
(see fig.l), and those that are constructed with three 3-phonon vertices oi Clo, (see fig. 2), 
that we call type II diagrams. The value of A is obtained after summing the contribution of 
all the Feynman diagrams depicted in fig. 1 and fig. 2, thus A = A^ + A^^. 

The scattering amplitude of diagram (i) of fig. 1 is given by 



iA] 



(i) 



3^3/2 



'^cliPa + Pb) ■ P{aPb} - {Pa + " 2u)p°p° - C^Pa ■ Pb) Gph {Pa + Pb^Pa + Pb) 

3c2(l - 2u) 



3c4 

-yiP^+Pb) -PidiPe-Pf}) + 2 

1 + 2u{5u - 4) - 2wp){pl + pMpyj] , 



{Pa + Pb)p\dPe ■ Pf} -{Pa+Pb)- P{dPePf} 

(4.3) 



while for diagram (i) of fig. 2 is 



3^3/2 



ctP 



2cl{pb - Pd) ■ p\bPd} - {Pb - - 2u)p°p° - c% ■ Pd) Gph [Pb - Pd^Pb - Pd) 

[cliPb - Pd) ■ {PaiPe + Pf) + (P° + p})Pa) " {pt " P^d) ((1 " 2n)p°(p° + pj) - C% ■ {Pe + Pf))] 
Qph {pI +P%Pe+ Pf) \2cl (pe + Pf) ■ p%Pf^ -{pl+p)){{l- 2u)pyf - C% ■ Pf) . 

(4.4) 
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Figure 1. Type I diagrams are formed with one 3-phonon vertex and one 4-phonon vertex of Clo- 



The curly brackets above indicate that the quantity has to be symmetrized with respect to 
the index inside the bracket, and is the phonon propagator 



(p0)2 



£2 

P 



(4.5) 



and at LO, we have Ep = CsP- All the remaining diagrams of fig. 1 and fig. 2 can be obtained 
from the expressions given above by relabeling the momenta and using crossing symmetry 
when necessary. 

Now it is important to realize that for certain configurations of the momenta the in- 
termediate propagators present in the diagrams on fig. 1 and fig. 2 will be on-shell, and 
if one considers that these are computed with a LO dispersion law, this would lead the 
corresponding amplitude to diverge. This fact has been previously recognized in Ref. [27]. 

It is possible to determine the configurations of the external momenta that lead the 
propagator to be on-shell by applying the kinematical considerations associated to the con- 
straints of energy and momentum conservation in every vertex. In all diagrams of both type 
I and type II the propagators are attached to a 3-phonon vertex with two external legs. 
When using the linear dispersion law it is easy to prove, see Appendix A, that the propaga- 
tor is on-shell when the momenta of the 3-phonon vertex are collinear. Therefore, for every 
diagram, the phase space contains regions where the propagators are on-shell. In type II di- 
agrams, since all verticies are of the 3-phonon type, all external momenta must be collinear. 
In type I diagrams, depending on the direction of the momentum flow in the propagator, in 
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Figure 2. Type II diagrams are constructed with three 3-phonon vertices of Clo 



the 4-phonon vertex we will have either one incoming and three outgoing phonons, or two 
incoming and two outgoing. When the propagator is on-shell, in the first case all external 
momenta will be collinear, see Appendix A, but not in the latter one. 

The existence of the mentioned collinear singularities in the computation of Fph is due 
to the fact that the LO Lagrangian used for its computation is not enough to describe the 
process under consideration. In order to cure these divergences, one should also consider the 
corrections associated to the NLO physics. We will not take into account NLO corrections 
to the different 3 and 4-phonon vertices, as following the power counting associated to the 
phonon EFT, it is easy to realize that these would only represent T^/A^ corrections to our 
results. When the phonon propagator is considered with a NLO dispersion law, and taking 
into account that 7 < 0, then energy and momentum conservation applied to every vertex 
of the different Feynman diagrams allows to deduce that there is no configuration of the 
external momenta that makes the internal phonon propagator to be on~shell (for details see 
Appendix A). 
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In the computation of Fph we will thus consider all the Feynman diagrams of fig. 1 
and fig. 2 using the phonon propagators with a NLO dispersion law. One then sees that, in 
this case, the almost collinear region (or small angle scattering region) of the available phase 
space is enhanced with respect to the rest of the phase space (or large angle scattering region) 
because the denominator of the phonon propagators in this region is small. In particular, 
this can be seen if we write the NLO phonon propagator as 

Sph {Pi +P°j,Pi+Pj) = i 

J4.6) 

Then, it is possible to see that in the almost collinear region, where the angle between pi and 
Pj is 9ij 0, the propagator behaves as ~ as compared to the region of large angle 

scattering, where the propagator behaves as ~ As a result, one can easily deduce that 

the almost collinear region of the phase space is enhanced with respect to the rest of the 
phase space. 

In our explicit numerical computations, and due to the fact that the NLO terms of 
the propagator are only relevant in the almost collinear region, we will approximate the 
propagator of eq. (4.6) as 

Gph {p- + P°j,Pi + Pj) ~ i [2clpiPj (1 - cos(%) - S^ipi + pjf)] . (4.7) 

Let us finally comment that some simple dimensional analysis allows us to see how the 
decay rate Fph associated to large or small angle scatterings behave [27]. We can obtain the 
dependence on T of the rate of the phonon number changing process by defining dimensionless 
variables associated to the momenta, Xi = Using these dimensionless variables, all the 
T dependence of the integral in eq. (4.1) factorizes. For the phase space d^^ the dependence 
is (see Appendix B). The T dependence of the amplitude corresponding to the diagrams 
in fig. 1 and fig. 2 can be computed by summing the number of legs of each one of the vertices 
of the diagram, which amounts to the total number of derivatives and hence the powers of 
the external momenta. To this number we have to substract the momentum dependence of 
the phonon propagator (s). The T dependence of the rate is then, two times the dependence 
of the amplitude added to the dependence of the phase space. For large angle collisions, then 
one sees that 

Fph cc . (4.8) 

In the almost collinear region the dependence of the propagators on the momentum is p'^ 
instead of p'^. Thus, the contributions to the rate coming from the collinear region have a 
different dependence on T. The T dependence of the collinear regions depends if we compute 
Fph with only the scattering matrix of the type I diagrams, or only with the scattering matrix 
of type II diagrams, or with the cross terms (that is, using ||^^^^^|| in Fph): 

r^ype I ^ 7^12 -ptypc II ,^8 pcross ^ rplO 

^ph o^-' ' -"-ph oc I , iph OCi 
4.1 Numerical results for the phonon decay rate 

The evaluation of the phonon rate depends on the features of the EoS as well as the value 
and density dependence of the neutron pairing gap in the core of neutron stars. While the 
APR nucleonic EoS is a common benchmark for all the EoS used in neutron star matter, 
it is still under debate the exact values of the ^Sq and neutron gaps and their density 



cliPi+Pj? 



1 + 27 



Pi + Pj 
Pi + Pj 



cl [Pi+Pj] 



l + 2-i{p,+pj) 



(4.9) 
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Figure 3. The ^Sq and averaged •^P2 neutron gaps as a function of the nucleon particle density in 
units of saturation density, uq — 0.16 fm~^. Two different neutron pairing gaps have been considered: 
a) ^So{A)+^P2{i) model, where the -^5*0 neutron gap is calculated in the BCS approach using different 
bare nucleon-nucleon interactions that converge towards a maximum neutron gap of about 3 MeV 
at pp sa 0.85fm~^ (parametrization A of Table I in Ref. [10]) while for we have taken the 
parametrization i (strong neutron superfluidity in the core); b) ^ So{a)+^ P2{h) model, where the 
^5*0 neutron gap incorporates medium polarization effects (parametrization a), whereas for the '^P2 
neutron gap we have considered the parametrization h (strong neutron superfluidity). 



Log r(cin ') 




Figure 4. Phonon rate in CGS units as a function of the temperature for different densities and the 
two neutron gap models. 

dependence [1]. It is, though, believed that ^5*0 neutron gap extends up to n ~ while 
neutrons are gapped in '^^2 well inside the core of the neutron star. 

The and averaged neutron gaps used throughout this work are shown in fig. 3. 
We have considered two very different gap models as a function of the density in order to illus- 
trate the model dependence of our results. Our first model, named hereafter ^So{A)+^P2{i), 
consists of the ^Sq neutron gap that results from the BCS approach using different bare 
nucleon-nucleon interactions that converge towards a maximum gap of about 3 MeV at 
Pf ~ 0.85fm~^ (parametrization A of Table I in Ref. [10]). The anisotropic '^^2 neutron 
gap is more challenging and not fully understood as one must extend BCS theory and cal- 
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Figure 5. The relative contribution of type I, type II and cross terms to the phonon rate as a function 
of the temperature for n = uq (left panel) and n = Auq (right panel). We use model ^So{A)+'^P2{i) for 
the neutron pairing gap. A similar behavior of these contributions is obtained for the ^So{a)+^P2{h) 
neutron gap model. 

culate several coupled equations while including relativistic effects since the gap extends for 
densities inside the core. We have taken the parametrization i (strong neutron superfluid- 
ity in the core) of Table I in Ref. [10] for the neutron angular averaged value, which 
presents a maximum value for the gap of approximately 1 MeV. The second model consid- 
ered, ^So{a)+^P2{h), goes beyond BCS for the ^Sq neutron gap as it incorporates medium 
polarization effects (parametrization a). The maximum value for the gap is then reduced to 
1 MeV. Moreover, for the neutron gap we have taken into account the parametrization 
h (strong neutron superfluidity) with a maximum value of about 0.5 MeV. We have, though, 
not considered weak neutron superfluidity in the core, as discussed in Ref. [10]. In this weak 
superfluid regime the values of the gap are A < 0.1 MeV for densities well inside the core. 
Thus, the corresponding transition temperatures from the superfluid to the normal phase 
are Tc ~ 1/2 A ^ 5 x 10^ K and no superfluidity is expected in the neutron star core for the 
temperatures studied. 

The phonon decay rate as a function of the temperature is displayed in fig. 4 for densities 
from 0.5no to 4no and for the two gap models previously discussed. These rates have been 
obtained solving eq. (4.1) numerically by means of Montecarlo integration. We can gauge 
the overall dependence of the rates with temperature by fitting a function of the form CT"^, 
where C and m are the parameters to be determined. The rates scale with T ~ 10^'^~^^K, 
the exact value depending on the density and the neutron gap model used. Given the 
previous dimensional analysis for the temperature dependence of the type I, type II and 
cross contributions to the rate, we conclude that the collinear region dominates over the 
non-collinear regime. 

Moreover, the relative importance of the type I, type II and cross contributions is shown 
in fig. 5 for n = uq and n = Auq using the ^50(^4) +^ -P2(0 neutron gap model. We observe 
that type II diagrams govern the behavior of the phonon decay rate up to T ~ lO^K from uq 
to 4no. As density increases, the dominance of type II terms extends to higher temperatures. 
Type I terms are sizeable as the temperature approaches to T ~ 10^'^~^'^K, specially for 
n = riQ, whereas the cross contributions remain small but non-negligible as temperature 
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Gap Model: 'So(A)+ 'P,(i) 
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Gap Model: '5o(a)+ T,© 
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Figure 6. Ci static bullc viscosity coefficient as a function of the temperature for various densities 
for the two neutron gap models. 
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Figure 7. (3 static bulk viscosity coefficient as a function of the temperature for various densities 
for the two neutron gap models. 

augments. A similar behaviom' is obtained for the type I, type II and cross contributions to 
the rate using the ^So{a) +^ P2{h) neutron gap modeh 



5 Values of the bulk viscosity coelRcients 

The static and frequency-dependent bulk viscosity coefficients of eqs. (3.4) and eqs. (3.9), 
respectively, result, on one hand, from the computation of the phonon rate and, on other 
hand, from the calculation of the Ii and I2 terms. The quantities /i,/2 depend on the EoS 
and, in particular, on the value and density dependence of the neutron pairing gaps (eqs. 3.7). 
Thus, the exact numerical results for the bulk viscosity coefficients calculated in the following 
will unavoidably depend on the chosen value of the neutron gaps and its density dependence, 
although the method of computation itself is rather general. 

In figs. 6,7 and 8 we present the static Cij Cs ^"^^ C2 bulk viscosity coefficients as a 
function of the temperature for densities between 0.5no and Auq. The bulk coefficients using 
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Figure 8. C2 static bulk viscosity coefficient as a function of the temperature for various densities 
for the two neutron gap models. 



the ^5o(^) P2{i) neutron gap model are shown on the left panels while those for the 
^5*0(0) +^ P2ih) neutron gap scheme are displayed on the right panels. As previously in- 
dicated, the C2 coefficient has the same physical meaning as the bulk viscosity in a normal 
fluid. We observe that the scaling with temperature of all three coefficients is given by the 
temperature dependence of the phonon decay rate. As we have seen in Sec. 4.1, the behavior 
with temperature of the rate is dominated by the contributions of type I, type II and cross 
diagrams stemming from the collinear region of the 2^-7-3 processes. In particular, type II 
terms govern the decay rate and, hence, the bulk viscosity coefficients up to T ~ lO^K for all 
the densities studied. On the other hand, each of the three coefficients has a distinct density 
dependence which results from the different combination of the density-dependent Ii and I2 
quantities, according to eqs. (3.4, 3.5). A slightly different behavior with density for the bulk 
coefficients is also manifest depending on the neutron gap model used. By comparing the 
results for both neutron gap schemes, we note that all bulk coefficients present bigger values 
for the ^S'o(a) +^ P2{h) case for all densities, partially due to the smaller phonon rates as 
seen in fig. 4. 

In fig. 9 we display the frequency-dependent ^2 coefficient at 4?io, which describes the 
damping of stellar pulsations with typical frequencies of w = 10^ — 10^s~^. When using the 
'^5o(^) +^ P2{i) neutron gap model, we observe on the left-hand side of fig. 9 that only for 
oj > 10^s~^ and T > lO^K we can distinguish between the static and frequency-dependent 
case. In fact, the left panel of fig. 10 shows that the static approximation is valid for the 
studied range of densities and temperatures, with the exception of 4no and T > lO^K since the 
characteristic frequency Uc becomes comparable to the typical value of the radial pulsations 
in stars. On the contrary, the ^2 coefficient is strongly dependent on the frequency if the 
^5*0(0) -|-^ P2{h) model is considered, as seen in the right-hand side of fig. 9. In this case the 
characteristic frequency of the right panel of fig. 10 becomes similar or even smaller than 
the typical stellar pulsation frequencies for all temperatures at n > 4no and, thus, the C2 
coefficient is suppressed with increasing frequency, an effect that can be easily inferred from 
eqs. (3.9,3.10). Similar behavior is expected for the frequency-dependent Ci and Ca bulk 
viscosity coefficients. 

We can now compare our results for the frequency-dependent C2 bulk viscosity coefficient 
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Figure 9. (2 frequency-dependent bulk viscosity coefficient as a function of the temperature for 
4no and frequencies between 10'^ — lO^s^^ for the two neutron gap models. 

coming from the collisions among superfluid phonons with the contribution stemming from 
direct URCA [8] and modified URCA [9] processes for a typical frequency of a; = 10^s~^. 
From the results presented in fig. 2 of Ref. [8] for the contributions of direct URCA processes 
at T = lO^K in non-superfluid matter, we see that the phonon contribution to the bulk 
coefficient using both neutron gap models is several orders of magnitude larger except in 
Model II [28] around densities of 2no, when the sudden opening of the URCA processes 
takes place. In fig. 7 the contributions to the viscosity from direct URCA processes in 
Model II at density 4no is plotted as a function of the temperature for non-superfluid and 
superfluid matter, showing that these contributions are smaller for all temperatures in the 
range 10^ < T{K) < 10^*^ although converging in size as the temperature increases. The 
comparison is similar for modified URCA processes. In fig. 1 of Ref. [9] these contributions 
in non-superfiuid matter as a function of the density for T = lO^K are displayed, being 
much smaller than the phonon contributions considering both neutron gap models except for 
densities about 2no in Model II when the URCA processes open up. In fig. 4 the viscosity 
for Model I [28] and 2no density is plotted as a function of the temperature for both non- 
superfiuid and superfluid matter, the contributions being much smaller at low temperatures 
and converging to the same order of magnitude at large temperatures. 

6 Conclusions 

We have computed the three bulk viscosity coefficients that appear in the superfiuid hy- 
drodynamic equations as arising from the collisions among phonons in superfluid neutron 
stars. We have presented a detailed analysis of how the phonon dispersion law determines 
the possible collisional processes relevant for the computation of these transport coefficients, 
and also how their explicit values depend on the EoS of the nucleonic matter inside the star 
as well as the neutron pairing gap. However, our method of computation is rather general, 
and could be used for different superfiuid systems, provided they share the same underly- 
ing symmetries. Only the knowledge of the EoS of the superfiuid and the specific form of 
the phonon dispersion law would be needed to extract the value of the three bulk viscosity 
coefficients from our general formulation. 
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Figure 10. Characteristic frequency uJc, defined in eq. (3.10), at different densities as a function of 
the temperature for the two neutron gap models. 



In this article we have used the APR EoS in a causal form [2] to describe the /3-stable 
nuclear matter inside the star, while two very distinct parameterizations of the ^Sq and the 
averaged neutron gaps have been considered [10]. Whereas the APR EoS is a common 
benchmark for all the nucleonic EoS, the neutron pairing gaps are still model dependent [1]. 
Thus, we have employed two very different neutron gap functions in order to test the model 
dependence of our results as the exact numerical values for the bulk viscosity coefficients 
depend unavoidably on this gap. Nevertheless, any future improvement in the determination 
of the neutron pairing gaps can be easily accommodated in our general scheme. 

Our results indicate that the dominant contribution of type I, type II and cross diagrams 
to the three bulk viscosity coefficients comes from the collinear regime of the 203 processes. 
In particular, the collinear behavior of the type II terms governs the temperature dependence 
of the bulk viscosity coefficients up to T ~ lO^K for all the densities studied. We have also 
analyzed the frequency-dependent bulk viscosity coefficients as compared to the static case. 
We find that it is possible to distinguish between static and frequency-dependent values for 
densities of n > 4no depending on the model used for the gap, being the frequency-dependent 
coefficient suppressed with respect to the static one as the frequency becomes bigger than 
the characteristic frequency for phonon collisions. Finally we have compared our results with 
those obtained for the bulk viscosities arising from direct and modified URCA processes [8, 9]. 
We conclude that, at T ~ lO^K and for typical radial pulsations of the star of w ~ 10^s~^, 
phonon collisions give the leading contribution to the bulk viscosities in the core, except 
for n ~ 2no when the sudden opening of the URCA processes take place. Note that our 
calculation for the bulk viscosity coefficients due to superfluid phonons is deemed to be valid 
up to Tc ~ 10^*^ K, which is the typical transition temperature to the normal fiuid. 

Our outcome can be used for studying damping of pulsations in neutron stars and grav- 
itational radiation driven instabilities in rotating neutron stars [29, 30], or the propagation 
of the sound waves within the star. In particular, it could be checked whether the phonon 
contribution to the bulk viscosities has any impact on the r-mode instability window of su- 
perfluid neutron stars, in the same way as it has been found that the phonon contribution 
to the shear viscosity modifies the r-mode instability window [16]. 
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A Allowed phonon collisions 

While using the phonon EFT at LO different phonon collisional processes are kinematically 
allowed, the physics beyond LO may introduce some kinematical restrictions on some scat- 
terings, as we discuss in this Appendix. Here we will restrict the discussion of the allowed 
kinematics when a phonon dispersion law of the form given in eq. (2.7) is considered. Then 
one can see that the sign of 7 plays a crucial role in determining whether some processes are 
or not possible. 

Let us first consider the decay of one to two phonons, or the reverse process, in general, 
1 -H- 2. Labelling the incoming particle as a and the outgoing as h and c, energy and 
momentum conservation imposes 

Ea = Eh + E,, (A.l) 

Pa = Pb+Pc- (A.2) 

Using the leading order dispersion relation, Ei = CsPi, in the energy conservation equation 
(A.l) and then using the momentum conservation relation (A.2) to eliminate pa from the 
former, we obtain an equality from which we can determine the value for the angle 9i,c between 
Pb and Pc 

PbPc{l-cos{ebc)) = 0. (A.3) 

The solutions pb = or pc = are naturally excluded, so one finds Obc = 0. Note that if we 
eliminate pb or pc instead of pa we will find Oac = ^afe = 0. We conclude that when all the 
legs of the 3-phonon vertex are on-shell the associated momenta are collinear. It is possible 
to compute the correction to 9bc due to the inclusion of NLO contributions to the dispersion 
relation. Lets define 69 as a small perturbation on the LO result, 9bc = + S9bc- Proceeding 
as we did before, we use the NLO dispersion relations in the energy conservation relation and 
eliminate pa using the momentum conservation. Then expanding both sides of the equation 
to first order in 7 and 69bc, we find the NLO correction to 9bc 

S9bc = V^{Pb + Pc) ■ (A.4) 

Obviously, for the one to two processes to be kinematically allowed, it is necessary that 7 > 0. 
In other words, when considering the NLO dispersion relations, 7 must be positive to have 
all the legs of the 3-phonon vertex on-shell. 

A similar analysis can be made for the one to three phonon process, 1^-7-3. Let us 
label the incoming particle as a and the outgoing particles as b, c and d. Solving momentum 
conservation equation for pa, and introducing it in the energy conservation equation while 
using the LO dispersion relation, we obtain 

PbPc{cos{9bc) - 1) + PbPd{cos{9bd) - 1) + PcPd{cos{9cd) - 1) = , (A. 5) 
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with 9ij being the angle between the vectors pi and pj. Excluding the solutions where two 
of the momenta are zero, we are left with the following solution 

Obc = 0, 0bd = O, dcd = 0. (A.6) 

If instead of solving for pa, we were to solve it for any of the outgoing momenta, this leads 
to 6ac = Oab = ^ad = 0. Now we cau study the corrections to the angles due to the NLO 
corrections. Following the same steps as in the 10 2 case, we define a small 69 correction 
to each one of the angles. We arrive to the following equation 

PbPcSelc + PbPdSO^d + PcPd^Ocd 

= 67 {Pb +Pc+ Pd) {PbPc + PbPd + plPb + plPd + PdPd + PdPc + 'iPbPcPd) , (A. 7) 

which only has solution for 7 > because all terms of the equation are positive. Thus the 
1^-7-3 process is kinematically allowed only when 7 > 0. This procedure for analysing the 
kinematics can be extended to any process of one to n phonons leading to the conclusion 
that these processes are only allowed for 7 > 0. 

To find a phonon number changing process that is kinematically allowed for 7 < we 
have to look at the two to three phonon processes, 2^3. We label the incoming particles 
a and b and the outgoing d, e, /. If we proceed as previously for the LO dispersion relation 
we obtain the following relation from the conservation laws 

PdPeil - COs{0de))+PdPfil - COs{0df)) + PeP/(l " COs(6'e/)) = 

PbPd{l - COs{6bd)) +PbPe{l - COs{9be)) - COs(6'fe/)) , 

which does not allow to determine any of the variables. Thus the two to three process is 
kinematically allowed regardless of the sign of 7. 

B Phase space integral 

We give here some details about the choice of variables in performing the phase space integral 
of eq. (4.1). We first integrate over d^Pb and pa making use of the momentum and the energy 
Dirac deltas respectively, which reduces the phase space to 

P*aPdPeP fdpfdpddpedQ.adQ.ddQ.edQ. f 

25(27r)llc6 {pd{l - COs{9ad)) +Pe{l " COs(0„e)) +Pf{l " COs(0,/))) ' ^ ' ' 

where dQ stand for the angular variables, and p* is defined as 

* _ PdPe(l - COs(6'rfe)) - COs{9df)) +PePf{l - COs{9ef)) 

Pd{l - COs{9ad)) + Pe(l - COs(6'ae)) + Pf{l " COs(6'a/)) 

We can further simplify the phase space expression by choosing a specific parametrization 
on the momenta. To define the orientation of a vector we need two angles. We have four 
vectors, so that amounts to eight angular variables. However, we have the freedom to choose 
the orientation of the reference frame. This freedom can be used to orientate the z axis along 
the direction of one of the momenta, and the zy plane to be parallel to the one generated by 
the same momentum and another of the momenta. We have chosen 

Pa = Pa (0, 0, 1) , 

Pd = Pdisin{9d)cos{(j)d), sin(6'd)sin(</)rf), cos{9d)) , 

Pe = Pe (0, Sin(6'e), COs(6'e)) , 

Pf = cos(0/), sm{9f)sm{(j)f), cos(0/)) . (B.3) 



(B.2) 
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Integrating the ciclic angular variables the phase space is reduced to 

p*aPdPeP f siniOd) sm{ee) sm{ef) M JJ, .IJ, ^ A A 

-ifi.o ^9 6 t — la w I 71 (a \\ \ n T^-^d9dd9edefd4>dd4>fdpddpedpf . 

16(27r)yco - cos(6'd)) +pe(l - cos(6le)) - cos(6'/))) 

(B.4) 
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